General project set-up

# Get all libraries and sources required to run the script
        library(plyr)
        library(dplyr)
        library(reshape2)

        library(ggplot2)
        library(ggthemes)
        library(scales)

        library(lme4)
        library(lmerTest)
        library(emmeans)
    
        library(skimr)

# Default ggplot settings

    Fill.colour<-scale_colour_manual(values = c("black", "gray70", "gray35"))

    ggthe_bw<-theme(plot.background=element_blank(),
          panel.grid.major.y = element_blank(),
          panel.grid.major.x = element_blank(),
          panel.grid.minor.x = element_blank(),
          panel.grid.minor.y = element_blank(),
          legend.box.background = element_rect(),
          panel.background =element_rect(fill = NA, color = "black")
          )+
    theme_bw()

1. DATA Exploration

1. Select qPCR info and define factors

# Import the data
  qPCR.variables  <- read.csv("Outputs/Ssid_SH_cell_ratio.csv", header = T)

# Remove blastate samples
  qPCR.blastate<-qPCR.variables[((qPCR.variables$Date=="2018-02-05") |
                                   (qPCR.variables$Date=="2018-03-08")), ]
  qPCR.variables <- droplevels(qPCR.variables[!rownames(qPCR.variables) %in% rownames(qPCR.blastate), ])  
  
 
# Variable types 
str(qPCR.variables)
## 'data.frame':    644 obs. of  17 variables:
##  $ Treatment : Factor w/ 3 levels "A","N","N+P": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Replicate : Factor w/ 2 levels "R1","R2": 1 1 1 1 1 1 1 1 1 1 ...
##  $ Fragment  : Factor w/ 162 levels "Ss_20-1","Ss_20-11",..: 27 30 31 48 49 52 55 58 97 100 ...
##  $ Genotype  : Factor w/ 7 levels "Ss_20","Ss_22",..: 2 2 2 2 3 3 3 3 5 5 ...
##  $ Date      : Factor w/ 5 levels "2017-11-15","2017-12-14",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Days      : int  0 0 0 0 0 0 0 0 0 0 ...
##  $ Time_Point: Factor w/ 5 levels "T13","T19","T22",..: 4 4 4 4 4 4 4 4 4 4 ...
##  $ Phase     : Factor w/ 3 levels "Baseline","Heat",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ Sample    : Factor w/ 644 levels "Ss_20-1_T5","Ss_20-1_T9",..: 109 124 128 188 191 204 217 232 380 393 ...
##  $ C.Ssid    : num  0 0 0 0 0 ...
##  $ D.Ssid    : num  0.0266 0.034 0.0239 0.0281 1.0265 ...
##  $ TotalSH   : num  0.0266 0.034 0.0239 0.0281 1.0265 ...
##  $ logC.SH   : num  NA NA NA NA NA ...
##  $ logD.SH   : num  -1.5756 -1.4689 -1.621 -1.5515 0.0114 ...
##  $ logSH     : num  -1.5756 -1.4689 -1.621 -1.5515 0.0114 ...
##  $ D.Prp     : num  1 1 1 1 1 ...
##  $ Community : Factor w/ 3 levels "C1","C3","D": 3 3 3 3 3 3 3 3 3 3 ...
  qPCR.variables$Genotype<-factor(qPCR.variables$Genotype, 
                                   levels=c("Ss_22","Ss_23","Ss_27", "Ss_28",
                                          "Ss_20", "Ss_24","Ss_30"))
  qPCR.variables$DaysF<-as.factor(qPCR.variables$Days)
  
  summary(qPCR.variables)
##  Treatment Replicate     Fragment    Genotype          Date    
##  A  :220   R1:333    Ss_20-17:  5   Ss_22:92   2017-11-15:140  
##  N  :228   R2:311    Ss_20-19:  5   Ss_23:96   2017-12-14:154  
##  N+P:196             Ss_20-22:  5   Ss_27:96   2018-02-01:151  
##                      Ss_20-26:  5   Ss_28:82   2018-02-23:103  
##                      Ss_20-29:  5   Ss_20:97   2018-03-06: 96  
##                      Ss_20-32:  5   Ss_24:85                   
##                      (Other) :614   Ss_30:96                   
##       Days        Time_Point       Phase              Sample   
##  Min.   :  0.00   T13:151    Baseline :140   Ss_20-1_T5  :  1  
##  1st Qu.: 29.00   T19:103    Heat     :199   Ss_20-1_T9  :  1  
##  Median : 78.00   T22: 96    Nutrients:305   Ss_20-11_T13:  1  
##  Mean   : 57.76   T5 :140                    Ss_20-11_T19:  1  
##  3rd Qu.:100.00   T9 :154                    Ss_20-11_T22:  1  
##  Max.   :111.00                              Ss_20-11_T9 :  1  
##                                              (Other)     :638  
##      C.Ssid             D.Ssid            TotalSH         
##  Min.   :0.000000   Min.   :0.000000   Min.   :0.0000034  
##  1st Qu.:0.000000   1st Qu.:0.000000   1st Qu.:0.0093391  
##  Median :0.002373   Median :0.006613   Median :0.0313298  
##  Mean   :0.046419   Mean   :0.059660   Mean   :0.1060796  
##  3rd Qu.:0.027225   3rd Qu.:0.044028   3rd Qu.:0.1286860  
##  Max.   :1.155796   Max.   :1.394495   Max.   :1.3944953  
##                                                           
##     logC.SH            logD.SH            logSH             D.Prp       
##  Min.   :-5.47264   Min.   :-4.7761   Min.   :-5.4726   Min.   :0.0000  
##  1st Qu.:-2.80684   1st Qu.:-2.2410   1st Qu.:-2.0297   1st Qu.:0.0000  
##  Median :-2.09749   Median :-1.6835   Median :-1.5040   Median :0.6792  
##  Mean   :-2.08904   Mean   :-1.7359   Mean   :-1.5221   Mean   :0.5242  
##  3rd Qu.:-1.27217   3rd Qu.:-1.0631   3rd Qu.:-0.8905   3rd Qu.:1.0000  
##  Max.   : 0.06288   Max.   : 0.1444   Max.   : 0.1444   Max.   :1.0000  
##  NA's   :184        NA's   :205                                         
##  Community DaysF    
##  C1:132    0  :140  
##  C3:171    29 :154  
##  D :341    78 :151  
##            100:103  
##            111: 96  
##                     
## 
  skim(qPCR.variables)
Data summary
Name qPCR.variables
Number of rows 644
Number of columns 18
_______________________
Column type frequency:
factor 10
numeric 8
________________________
Group variables None

Variable type: factor

skim_variable n_missing complete_rate ordered n_unique top_counts
Treatment 0 1 FALSE 3 N: 228, A: 220, N+P: 196
Replicate 0 1 FALSE 2 R1: 333, R2: 311
Fragment 0 1 FALSE 162 Ss_: 5, Ss_: 5, Ss_: 5, Ss_: 5
Genotype 0 1 FALSE 7 Ss_: 97, Ss_: 96, Ss_: 96, Ss_: 96
Date 0 1 FALSE 5 201: 154, 201: 151, 201: 140, 201: 103
Time_Point 0 1 FALSE 5 T9: 154, T13: 151, T5: 140, T19: 103
Phase 0 1 FALSE 3 Nut: 305, Hea: 199, Bas: 140
Sample 0 1 FALSE 644 Ss_: 1, Ss_: 1, Ss_: 1, Ss_: 1
Community 0 1 FALSE 3 D: 341, C3: 171, C1: 132
DaysF 0 1 FALSE 5 29: 154, 78: 151, 0: 140, 100: 103

Variable type: numeric

skim_variable n_missing complete_rate mean sd p0 p25 p50 p75 p100 hist
Days 0 1.00 57.76 41.59 0.00 29.00 78.00 100.00 111.00 ▆▆▁▆▇
C.Ssid 0 1.00 0.05 0.12 0.00 0.00 0.00 0.03 1.16 ▇▁▁▁▁
D.Ssid 0 1.00 0.06 0.14 0.00 0.00 0.01 0.04 1.39 ▇▁▁▁▁
TotalSH 0 1.00 0.11 0.17 0.00 0.01 0.03 0.13 1.39 ▇▁▁▁▁
logC.SH 184 0.71 -2.09 1.08 -5.47 -2.81 -2.10 -1.27 0.06 ▁▃▇▇▅
logD.SH 205 0.68 -1.74 0.94 -4.78 -2.24 -1.68 -1.06 0.14 ▁▂▆▇▃
logSH 0 1.00 -1.52 0.83 -5.47 -2.03 -1.50 -0.89 0.14 ▁▁▃▇▅
D.Prp 0 1.00 0.52 0.47 0.00 0.00 0.68 1.00 1.00 ▇▁▁▁▇

2. Exploratory graphs

  • Histograms
  HistoL_SH<-qplot(logSH, data=qPCR.variables, binwidth=0.15)
   HistoL_SH + facet_grid(Treatment~Date)

  ggplot(qPCR.variables, aes(logSH, fill = Treatment , colour = Treatment)) +
   geom_density(alpha = 0.1) + facet_wrap(~Date) + ggthe_bw

  • Log SH
logSHGenotype <- ggplot(qPCR.variables, aes(Days, logSH)) +
  geom_line(aes(colour=Fragment))+geom_point(aes(shape=factor(Community), colour=factor(Community)))+
  # geom_jitter(aes(colour=factor(Replicate))) +
      # stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2)+
      #stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5) +
      # stat_summary(fun.y=mean, geom="line", size =1, alpha=0.5) + 
      facet_grid(Treatment~Genotype) +
      ggthe_bw +theme(legend.position = "none" )

logSHGenotype + ylab("Relative log10 (S:H)") + xlab("Treatment") +  
      theme(axis.title.y=element_text(size=12), legend.position="none")

logSHTreatment <- ggplot(qPCR.variables, aes(Date, logSH, colour=factor(Treatment))) +
       # geom_jitter(aes(colour=factor(Treatment))) +
       stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2)+
       stat_summary(fun.y=mean, geom="line", size =1) +
       theme_bw() 
logSHTreatment + facet_wrap(~Genotype)

logSH_Replicate<- ggplot(qPCR.variables, aes (Days, logSH, 
                                             colour=factor(Replicate))) +
  stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2 )+
  stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5) +
  stat_summary(fun.y=mean, geom="line") + facet_grid (~Treatment) + ggthe_bw
logSH_Replicate 

logSH_Replicate+ facet_grid(~Community)

logSHTreatment <- ggplot (qPCR.variables, aes(Days, logSH, colour=Treatment)) + 
  ggtitle("A.") +  # geom_point(alpha=0.3) +
  ggthe_bw +  theme(legend.position="bottom") +
  
  stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", 
                    position = position_dodge(width = 5), width = 0.2 )+
  #stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5, 
       #            position = position_dodge(width = 1)) +
  stat_summary(fun.y=mean, geom="line", position = position_dodge(width = 5)) + 

  scale_y_continuous(breaks = seq(-5, 0.5, 0.5),
                    expand = c(0,0),
                  name=("log 10 (Total S/H cell ratio")) +
  scale_x_continuous(name="Days in the experiment",
                         limits = c(-2,113),
                         breaks = seq(0, 110, 15),  
                         expand = c(0, 0)) +  
 
      annotate("segment", x = 2, xend = 91, y = -3.5, yend = -3.5,
                colour = "gray90")+
      annotate("segment", x = 79, xend = 90, y = -3.5, yend = -3,
                colour = "gray90")+
      annotate("segment", x = 91, xend = 110, y = -3, yend = -3,
                colour = "gray90")
logSHTreatment

#logSHTreatment + facet_wrap(~Genotype)
logSHTreatment + facet_wrap(~Community)

  • Log C
logC_Treatment <- ggplot(qPCR.variables,
                         aes(Days, logC.SH, colour=factor(Fragment))) +
     geom_line()+ 
     #geom_jitter(aes(colour=factor(Replicate))) +
         #stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2)+
         #stat_summary(fun.y=mean, geom="line", size =1, alpha=0.5) + 
         ggthe_bw +  theme(legend.position = "none" )
logC_Treatment

logC_Treatment + facet_grid(Genotype~Treatment)

logC_Treatment + facet_grid(Community~Treatment)

  • Log D
logD_genotype <- ggplot(qPCR.variables,
                         aes(Days, logD.SH, colour=factor(Fragment))) +
     geom_line()+ 
     #geom_jitter(aes(colour=factor(Replicate))) +
         #stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2)+
         #stat_summary(fun.y=mean, geom="line", size =1, alpha=0.5) + 
         ggthe_bw +  theme(legend.position = "none" )
logD_genotype

logD_genotype + facet_grid(Genotype~Treatment)

logD_Treat <- ggplot(qPCR.variables,
                         aes(Days, logD.SH, colour=factor(Treatment))) +
     #geom_jitter(aes(colour=factor(Replicate))) +
        stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2)+
        stat_summary(fun.y=mean, geom="line", size =1, alpha=0.5) + 
        ggthe_bw +  theme(legend.position = "none" )
logD_Treat

logD_Treat + facet_grid(Community~Treatment)

*D proportion

qPCR.variables$D.Prp[qPCR.variables$Sample=="Ss_22-18_T22"]<-NA
qPCR.variables$Community[qPCR.variables$Sample=="Ss_22-18_T22"]<-"D"

qPCR.variables$D.Prp[qPCR.variables$Sample=="Ss_23-80_T19"]<-NA
qPCR.variables$Community[qPCR.variables$Sample=="Ss_23-80_T19"]<-"D"

D.PTreatment <- ggplot(qPCR.variables, aes(Days, D.Prp, colour=Fragment)) +
        ggthe_bw +  theme(legend.position="none") +  
        geom_jitter(aes(colour=factor(Fragment)), alpha=0.5, size=1) +
        geom_line() + facet_grid (Genotype~Treatment) +
  
      annotate("segment", x = 2, xend = 91, y = 0, yend = 0,
                colour = "gray90")+
  
      annotate("segment", x = 79, xend = 90, y = 0, yend = 0.3,
                colour = "gray90")+
      annotate("segment", x = 91, xend = 110, y = 0.3, yend = 0.3,
                colour = "gray90")+
  
        annotate("text", x = 8, y = 0.08, label = "(BL)", size=3, colour="gray")+
        annotate("text", x = 45, y = 0.08, label = "(C)", size=3, colour="gray")+
        annotate("text", x = 100, y = 0.08, label = "(H)", size=3, colour="gray")+  
  
      
      scale_y_continuous(breaks = seq(0, 1, 0.3),
                     name=("Durusdinium proportion (D/H)/(S/H)")) +
      scale_x_continuous(name="Days in the experiment",
                           limits = c(-2,113),
                         breaks = seq(0, 110, 45),  
                         expand = c(0, 0))
D.PTreatment

ggsave(file="Outputs/D.P_Treatment_Genotype.svg", plot=D.PTreatment, width=4.5, height=8)

2. Subset data

# Removing weird points and colonies not used

  qPCR.variables_2<-subset(qPCR.variables, Date!="2017-12-14")
  qPCR.variables_2<-subset(qPCR.variables_2, TotalSH<0.8)
  
  SH.0<-subset(qPCR.variables_2, Days<2)
  SH.C<-subset(qPCR.variables_2, Days<79)
    SH.C<-subset(SH.C, Days>77)
  SH.H<-subset(qPCR.variables_2, Days>80)
 
  
logSH_Genotype<- ggplot(qPCR.variables_2, aes (Days, logSH, 
                     colour=factor(Fragment))) +
  theme_bw() + geom_line()+ facet_grid (Genotype~Treatment)#+
logSH_Genotype + theme(legend.position = "none" ) + geom_point()

Log S/H graph

logSHTreatment <- ggplot (qPCR.variables_2, aes(Days, logSH, colour=Treatment)) +  theme_bw() + ggtitle("A.")+        Fill.colour+ ggthe_bw+
    theme(plot.background=element_blank(), 
            panel.border = element_blank(),
            axis.line = element_line(colour = "black"),
            panel.grid.major = element_blank(),
            panel.grid.minor = element_blank(),
            legend.position="bottom",
            strip.background = element_rect(fill="white")) +
       #geom_jitter(aes(colour=factor(Replicate))) +
       stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", 
                    position = position_dodge(width = 5), width = 0.2 )+
      #stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5, 
       #            position = position_dodge(width = 5)) +
       stat_summary(fun.y=mean, geom="line") + 
      
  annotate("segment", x = 2, xend = 91, y = -4.5, yend = -4.5,
                  colour = "gray90", linetype=1)+
        annotate("segment", x = 79, xend = 91, y = -4.5, yend = -4,
                  colour = "gray90", linetype=1)+
        annotate("segment", x = 91, xend = 110, y = -4, yend = -4,
                  colour = "gray90", linetype=1)+
        annotate("text", x = 45, y = -4.4, label = "Nutrients", size=3)+
        annotate("text", x = 99, y = -3.9, label = "Heat", size=3) +

      scale_y_continuous(breaks = seq(-5, 0.3, 0.5),
                    expand = c(0,0),
                   name=("Symbiodiniaceae / S. siderea (S/H) cell ratio")) +
      scale_x_continuous(name="Days in the experiment",
                         limits = c(-2,113),
                         breaks = seq(0, 110, 15),  
                         expand = c(0, 0))
  logSHTreatment 

  logSHTreatment + facet_grid(~Community)

  #ggsave(file="5.3A_AH.svg", plot=logSHTreatment, width=3.0, height=6)
  #ggsave(file="5.3A_AH_NoGenotype.svg", plot=logSHTreatment, width=3.0, height=3)

Log C/H graph

CHTreatment_2 <- ggplot (qPCR.variables_2, aes(Days, logC.SH, colour=Treatment)) +  theme_bw() + ggtitle("C.")+ Fill.colour+ ggthe_bw+
  annotate("segment", x = 2, xend = 91, y = -4.5, yend = -4.5,
                  colour = "gray90", linetype=1)+
        annotate("segment", x = 79, xend = 91, y = -4.5, yend = -4,
                  colour = "gray90", linetype=1)+
        annotate("segment", x = 91, xend = 110, y = -4, yend = -4,
                  colour = "gray90", linetype=1)+
        annotate("text", x = 45, y = -4.4, label = "Nutrients", size=3)+
        annotate("text", x = 99, y = -3.9, label = "Heat", size=3) +
       
       #geom_jitter(aes(colour=factor(Replicate))) +
       stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", 
                    position = position_dodge(width = 5), width = 0.2 )+
      #stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5, 
       #            position = position_dodge(width = 5)) +
       stat_summary(fun.y=mean, geom="line") +
      scale_x_continuous(name="Days in the experiment",
                         limits = c(-2,113),
                         breaks = seq(0, 110, 15),  
                         expand = c(0, 0))
  CHTreatment_2  

  #CHTreatment_2 + facet_grid (Genotype~Community)

  CHTreatment_2 + facet_grid(~Community)

  #ggsave(file="5.3A_AH_2.svg", plot=SHTreatment_2, width=3.0, height=6)
  #ggsave(file="5.3C_SH_NoGenotype_2.svg", plot=SHTreatment_2, width=3.0, height=3)

Log D/H graph

DHTreatment_2 <- ggplot (qPCR.variables_2, aes(Days, logD.SH, colour=Treatment)) +  theme_bw() + ggtitle("C.")+ Fill.colour+ ggthe_bw+
  annotate("segment", x = 2, xend = 91, y = -4.5, yend = -4.5,
                  colour = "gray90", linetype=1)+
        annotate("segment", x = 79, xend = 91, y = -4.5, yend = -4,
                  colour = "gray90", linetype=1)+
        annotate("segment", x = 91, xend = 110, y = -4, yend = -4,
                  colour = "gray90", linetype=1)+
        annotate("text", x = 45, y = -4.4, label = "Nutrients", size=3)+
        annotate("text", x = 99, y = -3.9, label = "Heat", size=3) +
       
       #geom_jitter(aes(colour=factor(Replicate))) +
       stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", 
                    position = position_dodge(width = 5), width = 0.2 )+
      #stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5, 
       #            position = position_dodge(width = 5)) +
       stat_summary(fun.y=mean, geom="line") +
      scale_x_continuous(name="Days in the experiment",
                         limits = c(-2,113),
                         breaks = seq(0, 110, 15),  
                         expand = c(0, 0))
  DHTreatment_2 + facet_grid (Genotype~Community)

  DHTreatment_2 + facet_grid(~Community)

D Proportion

D.PTreatment <- ggplot(qPCR.variables, aes(Days, D.Prp, colour=Fragment)) +
        ggthe_bw +  theme(legend.position="none") +  
        geom_jitter(aes(shape=factor(Treatment)), alpha=0.5) +
        geom_line()+ 
      annotate("segment", x = 2, xend = 91, y = 0, yend = 0,
                colour = "gray90")+
      #annotate("text", x = 45, y = 0.2, label = "Nutrients", size=3)+
    
      annotate("segment", x = 79, xend = 90, y = 0, yend = 0.3,
                colour = "gray90")+
    
      annotate("segment", x = 91, xend = 110, y = 0.3, yend = 0.3,
                colour = "gray90")+
      #annotate("text", x = 99, y = 0.3, label = "Heat", size=3)+

      scale_y_continuous(breaks = seq(0, 1, 0.3),
                     name=("Durusdinium proportion (D/H)/(S/H)")) +
      scale_x_continuous(name="Days in the experiment",
                           limits = c(-2,113),
                         breaks = seq(0, 110, 45),  
                         expand = c(0, 0))
      D.PTreatment + facet_grid (Genotype~Replicate)

D.PTreatment

3. GLM models

Baseline

LME_BaseLine<-lmer(logSH ~ Treatment + Community + (1|Genotype) + (1|Replicate), 
                     data=SH.0)
   step(LME_BaseLine)
## Backward reduced random-effect table:
## 
##                 Eliminated npar   logLik    AIC    LRT Df Pr(>Chisq)    
## <none>                        8  -75.016 166.03                         
## (1 | Genotype)           0    7  -82.009 178.02 13.986  1  0.0001842 ***
## (1 | Replicate)          0    7 -109.286 232.57 68.539  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##           Eliminated  Sum Sq Mean Sq NumDF   DenDF F value Pr(>F)
## Treatment          1 0.59627 0.29813     2 126.566  2.0641 0.1312
## Community          2 0.79889 0.39945     2   6.535  2.7209 0.1381
## 
## Model found:
## logSH ~ (1 | Genotype) + (1 | Replicate)
   anova(LME_BaseLine) # Treatments is not significant
   ranova(LME_BaseLine) # Treatments is significant!
   LME_BaseLine<-lmer(logSH ~ Genotype + (1|Replicate), 
                     data=SH.0)
   step(LME_BaseLine)
## Backward reduced random-effect table:
## 
##                 Eliminated npar   logLik    AIC    LRT Df Pr(>Chisq)    
## <none>                        9  -72.767 163.53                         
## (1 | Replicate)          0    8 -106.895 229.79 68.256  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##          Eliminated Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## Genotype          0 8.6374  1.4396     6   129  9.7555 7.505e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## logSH ~ Genotype + (1 | Replicate)
   anova(LME_BaseLine) # Treatments is not significant
   ranova(LME_BaseLine) # Treatments is significant!
# Multicomp
Ssid.SH.C.emm<-emmeans(LME_BaseLine, ~Genotype)
      #contrast(Ofav.SH.C.emm, "tukey")
    Ssid.SH.C_groups<-cld(Ssid.SH.C.emm, by=NULL) # compact-letter display
    Ssid.SH.C_groups
# Effect plot
  plot(emmeans(LME_BaseLine, ~Genotype), comparisons = TRUE) +
      coord_flip(xlim = NULL, ylim = NULL, expand = TRUE) + theme_bw()

logSHTreatmentBL <- ggplot (SH.0, aes(Genotype, logSH, colour=Genotype)) + 
  #ggtitle("A.") + 
  ggthe_bw +  theme(legend.position="bottom") +
  stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", 
                    position = position_dodge(width = 5), width = 0.2 )+
  stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5, 
                  position = position_dodge(width = 1)) +
  
  scale_y_continuous(breaks = seq(-2.5, 0, 0.5),
                    expand = c(0.1,0.1),
                  name=("log 10 (Total S/H cell ratio"))
    
logSHTreatmentBL

logSHTreatmentBL + facet_grid(~Community)

C: Effect of nutrient treatments at control temperature

LME_Ssid.C<-lmer(logSH ~ Treatment + Community + (1|Replicate) + (1|Genotype), 
                     data=SH.C)
     step (LME_Ssid.C) #  Remove Genotype
## Backward reduced random-effect table:
## 
##                 Eliminated npar   logLik    AIC    LRT Df Pr(>Chisq)    
## <none>                        8  -73.147 162.29                         
## (1 | Genotype)           1    7  -73.264 160.53  0.235  1     0.6278    
## (1 | Replicate)          0    6 -105.884 223.77 65.240  1   6.63e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##           Eliminated Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## Treatment          1 0.7873  0.3936     2   144  2.8627   0.06037 .  
## Community          0 9.4779  4.7390     2   146 33.6053 9.879e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## logSH ~ Community + (1 | Replicate)
     anova (LME_Ssid.C)
     ranova (LME_Ssid.C)
LME_Ssid.C<-lmer(logSH ~ Treatment + Community + (1|Replicate),
                     data=SH.C)
     step (LME_Ssid.C) #  Remove replicate
## Backward reduced random-effect table:
## 
##                 Eliminated npar   logLik    AIC   LRT Df Pr(>Chisq)    
## <none>                        7  -73.264 160.53                        
## (1 | Replicate)          0    6 -105.884 223.77 65.24  1   6.63e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##           Eliminated Sum Sq Mean Sq NumDF DenDF F value    Pr(>F)    
## Treatment          1 0.7873  0.3936     2   144  2.8627   0.06037 .  
## Community          0 9.4779  4.7390     2   146 33.6053 9.879e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## logSH ~ Community + (1 | Replicate)
     anova (LME_Ssid.C)
     ranova (LME_Ssid.C)
# Multicomp
Ssid.SH.C.emm<-emmeans(LME_Ssid.C, ~Community)
      #contrast(Ssid.SH.C.emm, "tukey")
    Ssid.SH.C_groups<-cld(Ssid.SH.C.emm, by=NULL) # compact-letter display
    Ssid.SH.C_groups
# Effect plot
  plot(emmeans(LME_Ssid.C, ~Treatment), comparisons = TRUE) +
      coord_flip(xlim = NULL, ylim = NULL, expand = TRUE) + theme_bw()

logSHTreatmentC <- ggplot (SH.C, aes(Treatment, logSH, colour=Treatment, shape=Treatment)) + 
  ggtitle("A.")+ Fill.colour +
  ggthe_bw +  theme(legend.position="bottom") +
  
  stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", 
                    position = position_dodge(width = 5), width = 0.2 )+
  stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5, 
                  position = position_dodge(width = 1)) +
  
  scale_y_continuous(breaks = seq(-2.5, 0, 0.5),
                    expand = c(0.1,0.1),
                  name=("log 10 (Total S/H cell ratio"))
    
logSHTreatmentC

logSHGenotypeC <- ggplot (SH.C, aes(Genotype, logSH, colour=Genotype)) + 
  ggthe_bw +  theme(legend.position="bottom") +
  
  stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", 
                    position = position_dodge(width = 5), width = 0.2 )+
  stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5, 
                  position = position_dodge(width = 1)) +
  
  scale_y_continuous(breaks = seq(-2.5, 0, 0.5),
                    expand = c(0.1,0.1),
                  name=("log 10 (C/H) cell ratio"))
    
logSHGenotypeC

H: Effect of pre-exposure to nutrient treatments during heat challenge

LME_Ssid.H<-lmer(logSH ~ Treatment * Community * DaysF+ (1|Replicate) + (1|Genotype), 
                     data=SH.H)
     step (LME_Ssid.H) #  Remove Replicate
## Backward reduced random-effect table:
## 
##                 Eliminated npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>                       21 -120.54 283.08                         
## (1 | Replicate)          1   20 -120.54 281.08  0.000  1          1    
## (1 | Genotype)           0   19 -157.51 353.02 73.949  1     <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                           Eliminated Sum Sq Mean Sq NumDF  DenDF F value
## Treatment:Community:DaysF          0 1.8668  0.4667     4 174.72  3.0217
##                            Pr(>F)  
## Treatment:Community:DaysF 0.01928 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## logSH ~ Treatment + Community + DaysF + (1 | Genotype) + Treatment:Community + 
##     Treatment:DaysF + Community:DaysF + Treatment:Community:DaysF
     anova (LME_Ssid.H)
     ranova (LME_Ssid.H)
LME_Ssid.H<-lmer(logSH ~ Treatment * Community * DaysF+ (1|Genotype),
                     data=SH.H)
     step (LME_Ssid.H) #  Remove replicate
## Backward reduced random-effect table:
## 
##                Eliminated npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>                      20 -120.54 281.08                         
## (1 | Genotype)          0   19 -157.51 353.02 73.949  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                           Eliminated Sum Sq Mean Sq NumDF  DenDF F value
## Treatment:Community:DaysF          0 1.8668  0.4667     4 174.72  3.0217
##                            Pr(>F)  
## Treatment:Community:DaysF 0.01928 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## logSH ~ Treatment * Community * DaysF + (1 | Genotype)
     anova (LME_Ssid.H)
     ranova (LME_Ssid.H)
# Multicomp
Ssid.SH.H.emm<-emmeans(LME_Ssid.H, ~Treatment*Community*DaysF)
      #contrast(Ssid.SH.H.emm, "tukey")
    Ssid.SH.H_groups<-cld(Ssid.SH.H.emm, by=NULL) # compact-letter display
    Ssid.SH.H_groups<-Ssid.SH.H_groups[order(Ssid.SH.H_groups$Day,
                           Ssid.SH.H_groups$Community, 
                           Ssid.SH.H_groups$Treatment), ]
    
    Ssid.SH.H_groups
# Effect plot
  plot(emmeans(LME_Ssid.H, ~DaysF*Community|Treatment), comparisons = TRUE) +
      coord_flip(xlim = NULL, ylim = NULL, expand = TRUE) + theme_bw()

logSHTreatmentH <- ggplot (SH.H, aes(Treatment, logSH, colour=Treatment, shape=Treatment)) + 
  ggtitle("A.")+ Fill.colour +
  ggthe_bw +  theme(legend.position="bottom") +
  
  stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", 
                    position = position_dodge(width = 5), width = 0.2 )+
  stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5, 
                  position = position_dodge(width = 1)) +
  
  scale_y_continuous(breaks = seq(-3, 0, 0.5),
                    expand = c(0.1,0.1),
                  name=("log 10 (Total S/H cell ratio"))
    
logSHTreatmentH

logSHTreatmentH + facet_grid(~Community)

logSHTreatmentH + facet_grid(DaysF~Community) #+  geom_jitter()

logSHTreatmentH <- ggplot (SH.H, aes(DaysF, logSH)) + 
  ggtitle("A.")+ Fill.colour +
  ggthe_bw +  theme(legend.position="bottom") +
  
  stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", 
                    position = position_dodge(width = 5), width = 0.2 )+
  stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5, 
                  position = position_dodge(width = 1)) +
  
  scale_y_continuous(breaks = seq(-3, 0, 0.5),
                    expand = c(0.1,0.1),
                  name=("log 10 (Total S/H cell ratio"))
    
logSHTreatmentH

logSHTreatmentH + facet_grid(~Community)

logSHTreatmentH + facet_grid(Treatment~Community) #+  geom_jitter()

logSHGenotypeH <- ggplot (SH.H, aes(Community, logSH, colour=Genotype)) + 
  ggthe_bw +  theme(legend.position="bottom") +
  
  stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 0.2 )+
  stat_summary(fun.y=mean, geom="point", size =3, alpha=0.5) +
  
  scale_y_continuous(breaks = seq(-2.5, 0, 0.5),
                    expand = c(0.1,0.1),
                  name=("log 10 (Total S/H cell ratio"))
    
logSHGenotypeH

logSHGenotypeH + facet_grid(~DaysF)

All phases

  • Dasys as a factor
LM_Ssid.qpCR<-lmer(logSH ~ Treatment * DaysF * Community + (1|Replicate) + (1|Genotype/Fragment), 
                   data= qPCR.variables_2)
    step(LM_Ssid.qpCR)
## Backward reduced random-effect table:
## 
##                         Eliminated npar  logLik    AIC    LRT Df
## <none>                               40 -347.95 775.90          
## (1 | Replicate)                  1   39 -347.95 773.90  0.000  1
## (1 | Fragment:Genotype)          2   38 -347.95 771.90  0.000  1
## (1 | Genotype)                   0   37 -377.52 829.03 59.127  1
##                         Pr(>Chisq)    
## <none>                                
## (1 | Replicate)                  1    
## (1 | Fragment:Genotype)          1    
## (1 | Genotype)           1.478e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                           Eliminated  Sum Sq Mean Sq NumDF  DenDF F value
## Treatment:DaysF:Community          1  2.0375  0.1698    12 443.32  0.7963
## Treatment:Community                2  0.8016  0.2004     4 455.42  0.9449
## Treatment:DaysF                    0  5.1227  0.8538     6 459.60  4.0277
## DaysF:Community                    0 24.7264  4.1211     6 459.71 19.4411
##                              Pr(>F)    
## Treatment:DaysF:Community 0.6545075    
## Treatment:Community       0.4376783    
## Treatment:DaysF           0.0006045 ***
## DaysF:Community           < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## logSH ~ Treatment + DaysF + Community + (1 | Genotype) + Treatment:DaysF + 
##     DaysF:Community
    anova(LM_Ssid.qpCR)
    ranova(LM_Ssid.qpCR)
LM_Ssid.qpCR<-lmer(logSH ~ Treatment + DaysF + Community + 
                     (1 | Genotype) + 
                     Treatment:DaysF +  DaysF:Community,
                   data= qPCR.variables_2)
      step(LM_Ssid.qpCR)
## Backward reduced random-effect table:
## 
##                Eliminated npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>                      22 -346.54 737.08                         
## (1 | Genotype)          0   21 -376.97 795.94 60.855  1  6.144e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                 Eliminated  Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)
## Treatment:DaysF          0  5.1227  0.8538     6 459.60  4.0277 0.0006045
## DaysF:Community          0 24.7264  4.1211     6 459.71 19.4411 < 2.2e-16
##                    
## Treatment:DaysF ***
## DaysF:Community ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## logSH ~ Treatment + DaysF + Community + (1 | Genotype) + Treatment:DaysF + 
##     DaysF:Community
      anova(LM_Ssid.qpCR)
      ranova(LM_Ssid.qpCR)
      summary(LM_Ssid.qpCR)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## logSH ~ Treatment + DaysF + Community + (1 | Genotype) + Treatment:DaysF +  
##     DaysF:Community
##    Data: qPCR.variables_2
## 
## REML criterion at convergence: 693.1
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3953 -0.6798 -0.0318  0.7153  3.3216 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Genotype (Intercept) 0.4236   0.6508  
##  Residual             0.2120   0.4604  
## Number of obs: 486, groups:  Genotype, 7
## 
## Fixed effects:
##                        Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)            -1.89344    0.28407   8.84781  -6.665 9.97e-05 ***
## TreatmentN             -0.07281    0.09363 459.23851  -0.778 0.437183    
## TreatmentN+P           -0.23487    0.09896 459.27507  -2.373 0.018035 *  
## DaysF78                 0.18335    0.14590 459.31051   1.257 0.209516    
## DaysF100               -0.89065    0.15375 459.46923  -5.793 1.28e-08 ***
## DaysF111               -1.12847    0.15765 459.46253  -7.158 3.27e-12 ***
## CommunityC3             1.42388    0.22411 387.61081   6.353 5.91e-10 ***
## CommunityD              0.14326    0.14828 465.63001   0.966 0.334479    
## TreatmentN:DaysF78      0.09799    0.13142 459.23892   0.746 0.456285    
## TreatmentN+P:DaysF78    0.08542    0.13575 459.26579   0.629 0.529517    
## TreatmentN:DaysF100     0.48243    0.14265 459.33134   3.382 0.000781 ***
## TreatmentN+P:DaysF100   0.44027    0.15202 459.29481   2.896 0.003957 ** 
## TreatmentN:DaysF111     0.09349    0.14708 460.66389   0.636 0.525338    
## TreatmentN+P:DaysF111   0.46388    0.15488 459.43015   2.995 0.002891 ** 
## DaysF78:CommunityC3     0.35719    0.16195 459.31568   2.206 0.027907 *  
## DaysF100:CommunityC3   -0.11623    0.17443 459.35517  -0.666 0.505507    
## DaysF111:CommunityC3   -1.21190    0.19904 459.31114  -6.089 2.40e-09 ***
## DaysF78:CommunityD      0.10888    0.14638 459.32666   0.744 0.457370    
## DaysF100:CommunityD     0.03883    0.15790 459.35328   0.246 0.805831    
## DaysF111:CommunityD     0.26168    0.15986 460.69394   1.637 0.102321    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Emmeans
Ssid.YII.emm<-emmeans(LM_Ssid.qpCR, ~Treatment*DaysF*Community)
     #contrast(Ssid.YII.emm, "tukey")

# Pairwise comparisons
    #pairs(Ssid.YII.emm) # same than contrast(Ssid.YII.emm, "tukey")
  plot(emmeans(LM_Ssid.qpCR, ~Treatment|DaysF), comparisons = TRUE) +
        theme_bw() # Tukey comparission (do not trust CI to compare EMMs)

    plot(emmeans(LM_Ssid.qpCR, ~Treatment|DaysF|Community), comparisons = TRUE) +
      coord_flip(xlim = NULL, ylim = NULL, expand = TRUE) + facet_grid(Community~DaysF) +
      theme_bw()

  #CLD  
  SH_Multicomp<-cld(Ssid.YII.emm, by=NULL) # compact-letter display
      SH_Multicomp<-SH_Multicomp[order(SH_Multicomp$DaysF, 
                                       SH_Multicomp$Community,
                                       SH_Multicomp$Treatment), ]
    #write.csv(SH_Multicomp, "Outputs/Ssid_SH_Multicomp.csv")  
 # 2. Predict values:
    pred_Ssid1 <- predict(LM_Ssid.qpCR, re.form = NA)

  #3. Bootstrap CI:
    Ss.boot1 <- bootMer(LM_Ssid.qpCR, predict, nsim = 1000, re.form = NULL) # include random effects, reduce CI lot!
    std.err <- apply(Ss.boot1$t, 2, sd)
    CI.lo_1 <- pred_Ssid1 - std.err*1.96
    CI.hi_1 <- pred_Ssid1 + std.err*1.96

  #Plot
  Model_Ss_1b_plot<- ggplot(
    qPCR.variables_2, aes(x = Days, y = logSH, colour = Treatment)) +
    geom_line(aes(y = pred_Ssid1),size=2) +
    #geom_point(aes(fill=factor(Treatment)),
    #         shape = 21, colour = "black", size = 2, stroke = 0.3, alpha=0.3) +
    geom_ribbon(aes(ymin = CI.lo_1, ymax = CI.hi_1),
                size=2, alpha = 0.1, linetype = 0) +
    #scale_color_manual(values=my_colours) +
    #scale_fill_manual(values=my_colours) +
    #scale_y_continuous(name=expression(~italic("Fv / Fm")),
    #                   limits = c(0.5, 0.61), 
    #                   breaks = seq(0.4, 0.6, by=0.1), expand = c(0,0))+
    #scale_x_continuous("Days in the experiment", limits = c(0, 78),
    #                 breaks = seq(0, 76, by=7), expand = c(0,0))+
    
     stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1,
                   position = position_dodge(1) )+
    stat_summary(fun.y=mean, geom="line", position = position_dodge(1), 
                  linetype=1, alpha=1)
    # stat_summary(fun.y=mean, geom="point", size =1,
    #             position=position_dodge(width=1), alpha=0.8)  +
    ggthe_bw
## List of 69
##  $ legend.box.background     : list()
##   ..- attr(*, "class")= chr [1:2] "element_blank" "element"
##  $ panel.background          :List of 5
##   ..$ fill         : chr "white"
##   ..$ colour       : logi NA
##   ..$ size         : NULL
##   ..$ linetype     : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_rect" "element"
##  $ panel.grid.major.x        : list()
##   ..- attr(*, "class")= chr [1:2] "element_blank" "element"
##  $ panel.grid.major.y        : list()
##   ..- attr(*, "class")= chr [1:2] "element_blank" "element"
##  $ panel.grid.minor.x        : list()
##   ..- attr(*, "class")= chr [1:2] "element_blank" "element"
##  $ panel.grid.minor.y        : list()
##   ..- attr(*, "class")= chr [1:2] "element_blank" "element"
##  $ plot.background           :List of 5
##   ..$ fill         : NULL
##   ..$ colour       : chr "white"
##   ..$ size         : NULL
##   ..$ linetype     : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_rect" "element"
##  $ line                      :List of 6
##   ..$ colour       : chr "black"
##   ..$ size         : num 0.5
##   ..$ linetype     : num 1
##   ..$ lineend      : chr "butt"
##   ..$ arrow        : logi FALSE
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_line" "element"
##  $ rect                      :List of 5
##   ..$ fill         : chr "white"
##   ..$ colour       : chr "black"
##   ..$ size         : num 0.5
##   ..$ linetype     : num 1
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_rect" "element"
##  $ text                      :List of 11
##   ..$ family       : chr ""
##   ..$ face         : chr "plain"
##   ..$ colour       : chr "black"
##   ..$ size         : num 11
##   ..$ hjust        : num 0.5
##   ..$ vjust        : num 0.5
##   ..$ angle        : num 0
##   ..$ lineheight   : num 0.9
##   ..$ margin       : 'margin' num [1:4] 0pt 0pt 0pt 0pt
##   .. ..- attr(*, "valid.unit")= int 8
##   .. ..- attr(*, "unit")= chr "pt"
##   ..$ debug        : logi FALSE
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.title.x              :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : NULL
##   ..$ vjust        : num 1
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 2.75pt 0pt 0pt 0pt
##   .. ..- attr(*, "valid.unit")= int 8
##   .. ..- attr(*, "unit")= chr "pt"
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.title.x.top          :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : NULL
##   ..$ vjust        : num 0
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 0pt 0pt 2.75pt 0pt
##   .. ..- attr(*, "valid.unit")= int 8
##   .. ..- attr(*, "unit")= chr "pt"
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.title.y              :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : NULL
##   ..$ vjust        : num 1
##   ..$ angle        : num 90
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 0pt 2.75pt 0pt 0pt
##   .. ..- attr(*, "valid.unit")= int 8
##   .. ..- attr(*, "unit")= chr "pt"
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.title.y.right        :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : NULL
##   ..$ vjust        : num 0
##   ..$ angle        : num -90
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 0pt 0pt 0pt 2.75pt
##   .. ..- attr(*, "valid.unit")= int 8
##   .. ..- attr(*, "unit")= chr "pt"
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.text                 :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : chr "grey30"
##   ..$ size         : 'rel' num 0.8
##   ..$ hjust        : NULL
##   ..$ vjust        : NULL
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : NULL
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.text.x               :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : NULL
##   ..$ vjust        : num 1
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 2.2pt 0pt 0pt 0pt
##   .. ..- attr(*, "valid.unit")= int 8
##   .. ..- attr(*, "unit")= chr "pt"
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.text.x.top           :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : NULL
##   ..$ vjust        : num 0
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 0pt 0pt 2.2pt 0pt
##   .. ..- attr(*, "valid.unit")= int 8
##   .. ..- attr(*, "unit")= chr "pt"
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.text.y               :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : num 1
##   ..$ vjust        : NULL
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 0pt 2.2pt 0pt 0pt
##   .. ..- attr(*, "valid.unit")= int 8
##   .. ..- attr(*, "unit")= chr "pt"
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.text.y.right         :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : num 0
##   ..$ vjust        : NULL
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 0pt 0pt 0pt 2.2pt
##   .. ..- attr(*, "valid.unit")= int 8
##   .. ..- attr(*, "unit")= chr "pt"
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ axis.ticks                :List of 6
##   ..$ colour       : chr "grey20"
##   ..$ size         : NULL
##   ..$ linetype     : NULL
##   ..$ lineend      : NULL
##   ..$ arrow        : logi FALSE
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_line" "element"
##  $ axis.ticks.length         : 'unit' num 2.75pt
##   ..- attr(*, "valid.unit")= int 8
##   ..- attr(*, "unit")= chr "pt"
##  $ axis.ticks.length.x       : NULL
##  $ axis.ticks.length.x.top   : NULL
##  $ axis.ticks.length.x.bottom: NULL
##  $ axis.ticks.length.y       : NULL
##  $ axis.ticks.length.y.left  : NULL
##  $ axis.ticks.length.y.right : NULL
##  $ axis.line                 : list()
##   ..- attr(*, "class")= chr [1:2] "element_blank" "element"
##  $ axis.line.x               : NULL
##  $ axis.line.y               : NULL
##  $ legend.background         :List of 5
##   ..$ fill         : NULL
##   ..$ colour       : logi NA
##   ..$ size         : NULL
##   ..$ linetype     : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_rect" "element"
##  $ legend.margin             : 'margin' num [1:4] 5.5pt 5.5pt 5.5pt 5.5pt
##   ..- attr(*, "valid.unit")= int 8
##   ..- attr(*, "unit")= chr "pt"
##  $ legend.spacing            : 'unit' num 11pt
##   ..- attr(*, "valid.unit")= int 8
##   ..- attr(*, "unit")= chr "pt"
##  $ legend.spacing.x          : NULL
##  $ legend.spacing.y          : NULL
##  $ legend.key                :List of 5
##   ..$ fill         : chr "white"
##   ..$ colour       : logi NA
##   ..$ size         : NULL
##   ..$ linetype     : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_rect" "element"
##  $ legend.key.size           : 'unit' num 1.2lines
##   ..- attr(*, "valid.unit")= int 3
##   ..- attr(*, "unit")= chr "lines"
##  $ legend.key.height         : NULL
##  $ legend.key.width          : NULL
##  $ legend.text               :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : 'rel' num 0.8
##   ..$ hjust        : NULL
##   ..$ vjust        : NULL
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : NULL
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ legend.text.align         : NULL
##  $ legend.title              :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : num 0
##   ..$ vjust        : NULL
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : NULL
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ legend.title.align        : NULL
##  $ legend.position           : chr "right"
##  $ legend.direction          : NULL
##  $ legend.justification      : chr "center"
##  $ legend.box                : NULL
##  $ legend.box.margin         : 'margin' num [1:4] 0cm 0cm 0cm 0cm
##   ..- attr(*, "valid.unit")= int 1
##   ..- attr(*, "unit")= chr "cm"
##  $ legend.box.spacing        : 'unit' num 11pt
##   ..- attr(*, "valid.unit")= int 8
##   ..- attr(*, "unit")= chr "pt"
##  $ panel.border              :List of 5
##   ..$ fill         : logi NA
##   ..$ colour       : chr "grey20"
##   ..$ size         : NULL
##   ..$ linetype     : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_rect" "element"
##  $ panel.spacing             : 'unit' num 5.5pt
##   ..- attr(*, "valid.unit")= int 8
##   ..- attr(*, "unit")= chr "pt"
##  $ panel.spacing.x           : NULL
##  $ panel.spacing.y           : NULL
##  $ panel.grid                :List of 6
##   ..$ colour       : chr "grey92"
##   ..$ size         : NULL
##   ..$ linetype     : NULL
##   ..$ lineend      : NULL
##   ..$ arrow        : logi FALSE
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_line" "element"
##  $ panel.grid.minor          :List of 6
##   ..$ colour       : NULL
##   ..$ size         : 'rel' num 0.5
##   ..$ linetype     : NULL
##   ..$ lineend      : NULL
##   ..$ arrow        : logi FALSE
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_line" "element"
##  $ panel.ontop               : logi FALSE
##  $ plot.title                :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : 'rel' num 1.2
##   ..$ hjust        : num 0
##   ..$ vjust        : num 1
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 0pt 0pt 5.5pt 0pt
##   .. ..- attr(*, "valid.unit")= int 8
##   .. ..- attr(*, "unit")= chr "pt"
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ plot.subtitle             :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : num 0
##   ..$ vjust        : num 1
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 0pt 0pt 5.5pt 0pt
##   .. ..- attr(*, "valid.unit")= int 8
##   .. ..- attr(*, "unit")= chr "pt"
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ plot.caption              :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : 'rel' num 0.8
##   ..$ hjust        : num 1
##   ..$ vjust        : num 1
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 5.5pt 0pt 0pt 0pt
##   .. ..- attr(*, "valid.unit")= int 8
##   .. ..- attr(*, "unit")= chr "pt"
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ plot.tag                  :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : 'rel' num 1.2
##   ..$ hjust        : num 0.5
##   ..$ vjust        : num 0.5
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : NULL
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ plot.tag.position         : chr "topleft"
##  $ plot.margin               : 'margin' num [1:4] 5.5pt 5.5pt 5.5pt 5.5pt
##   ..- attr(*, "valid.unit")= int 8
##   ..- attr(*, "unit")= chr "pt"
##  $ strip.background          :List of 5
##   ..$ fill         : chr "grey85"
##   ..$ colour       : chr "grey20"
##   ..$ size         : NULL
##   ..$ linetype     : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_rect" "element"
##  $ strip.placement           : chr "inside"
##  $ strip.text                :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : chr "grey10"
##   ..$ size         : 'rel' num 0.8
##   ..$ hjust        : NULL
##   ..$ vjust        : NULL
##   ..$ angle        : NULL
##   ..$ lineheight   : NULL
##   ..$ margin       : 'margin' num [1:4] 4.4pt 4.4pt 4.4pt 4.4pt
##   .. ..- attr(*, "valid.unit")= int 8
##   .. ..- attr(*, "unit")= chr "pt"
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ strip.text.x              : NULL
##  $ strip.text.y              :List of 11
##   ..$ family       : NULL
##   ..$ face         : NULL
##   ..$ colour       : NULL
##   ..$ size         : NULL
##   ..$ hjust        : NULL
##   ..$ vjust        : NULL
##   ..$ angle        : num -90
##   ..$ lineheight   : NULL
##   ..$ margin       : NULL
##   ..$ debug        : NULL
##   ..$ inherit.blank: logi TRUE
##   ..- attr(*, "class")= chr [1:2] "element_text" "element"
##  $ strip.switch.pad.grid     : 'unit' num 2.75pt
##   ..- attr(*, "valid.unit")= int 8
##   ..- attr(*, "unit")= chr "pt"
##  $ strip.switch.pad.wrap     : 'unit' num 2.75pt
##   ..- attr(*, "valid.unit")= int 8
##   ..- attr(*, "unit")= chr "pt"
##  - attr(*, "class")= chr [1:2] "theme" "gg"
##  - attr(*, "complete")= logi TRUE
##  - attr(*, "validate")= logi TRUE
  Model_Ss_1b_plot + facet_wrap(~Community)

  • Days as continous
LM_Ssid.qpCR2<-lmer(logSH ~ Treatment * Days * Community + (1|Replicate) + (1|Genotype/Fragment), 
                   data= qPCR.variables_2)
    step(LM_Ssid.qpCR2)
## Backward reduced random-effect table:
## 
##                         Eliminated npar  logLik    AIC   LRT Df Pr(>Chisq)
## <none>                               22 -568.24 1180.5                    
## (1 | Replicate)                  1   21 -568.24 1178.5  0.00  1          1
## (1 | Fragment:Genotype)          2   20 -568.24 1176.5  0.00  1          1
## (1 | Genotype)                   0   19 -589.20 1216.4 41.93  1  9.461e-11
##                            
## <none>                     
## (1 | Replicate)            
## (1 | Fragment:Genotype)    
## (1 | Genotype)          ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                          Eliminated  Sum Sq Mean Sq NumDF  DenDF F value
## Treatment:Days:Community          1  1.7016  0.4254     4 461.15  0.9054
## Days:Community                    2  1.7500  0.8750     2 465.80  1.8642
## Treatment:Community               3  3.2309  0.8077     4 467.27  1.7133
## Community                         0 22.9355 11.4678     2 364.82 24.1860
## Treatment:Days                    0  3.0544  1.5272     2 471.54  3.2209
##                             Pr(>F)    
## Treatment:Days:Community    0.4606    
## Days:Community              0.1562    
## Treatment:Community         0.1458    
## Community                1.369e-10 ***
## Treatment:Days              0.0408 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## logSH ~ Treatment + Days + Community + (1 | Genotype) + Treatment:Days
    anova(LM_Ssid.qpCR2)
    ranova(LM_Ssid.qpCR2)
LM_Ssid.qpCR2<-lmer(logSH ~Treatment + Days + Community + (1|Genotype) + Treatment:Days,
                   data= qPCR.variables_2)
      step(LM_Ssid.qpCR2)
## Backward reduced random-effect table:
## 
##                Eliminated npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>                      10 -542.81 1105.6                         
## (1 | Genotype)          0    9 -564.04 1146.1 42.471  1  7.173e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                Eliminated  Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)
## Community               0 22.9355 11.4678     2 364.82 24.1860 1.369e-10
## Treatment:Days          0  3.0544  1.5272     2 471.54  3.2209    0.0408
##                   
## Community      ***
## Treatment:Days *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## logSH ~ Treatment + Days + Community + (1 | Genotype) + Treatment:Days
      anova(LM_Ssid.qpCR2)
      ranova(LM_Ssid.qpCR2)
      summary(LM_Ssid.qpCR2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## logSH ~ Treatment + Days + Community + (1 | Genotype) + Treatment:Days
##    Data: qPCR.variables_2
## 
## REML criterion at convergence: 1085.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.5286 -0.6117 -0.1087  0.6507  2.3774 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Genotype (Intercept) 0.6803   0.8248  
##  Residual             0.4741   0.6886  
## Number of obs: 486, groups:  Genotype, 7
## 
## Fixed effects:
##                     Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)        -1.853179   0.360321   8.344461  -5.143 0.000773 ***
## TreatmentN         -0.096926   0.136250 471.094970  -0.711 0.477201    
## TreatmentN+P       -0.266358   0.143608 471.050405  -1.855 0.064257 .  
## Days               -0.008168   0.001209 471.227691  -6.757 4.18e-11 ***
## CommunityC3         1.818818   0.283439 290.421477   6.417 5.64e-10 ***
## CommunityD          0.262047   0.169203 453.518949   1.549 0.122148    
## TreatmentN:Days     0.003229   0.001701 471.924639   1.898 0.058301 .  
## TreatmentN+P:Days   0.004293   0.001799 471.108284   2.386 0.017420 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmnN TrtN+P Days   CmmnC3 CmmntD TrtN:D
## TreatmentN  -0.183                                          
## TreatmntN+P -0.185  0.480                                   
## Days        -0.248  0.591  0.564                            
## CommunityC3 -0.371 -0.026  0.005  0.053                     
## CommunityD  -0.382 -0.009  0.011  0.047  0.581              
## TrtmntN:Dys  0.149 -0.834 -0.400 -0.706  0.063 -0.018       
## TrtmntN+P:D  0.141 -0.399 -0.839 -0.669  0.021  0.020  0.477
 # 2. Predict values:
    pred_Ssid1 <- predict(LM_Ssid.qpCR2,re.form = NA)

  #3. Bootstrap CI:
    Ss.boot1 <- bootMer(LM_Ssid.qpCR2, predict, nsim = 1000, re.form = NULL) # include random effects, reduce CI lot!
    std.err <- apply(Ss.boot1$t, 2, sd)
    CI.lo_1 <- pred_Ssid1 - std.err*1.96
    CI.hi_1 <- pred_Ssid1 + std.err*1.96

  #Plot
  Model_Ss_1b_plot<- ggplot(
    qPCR.variables_2, aes(x = Days, y = logSH, colour = Treatment)) +
    geom_line(aes(y = pred_Ssid1),size=2) +
    #geom_point(aes(fill=factor(Treatment)),
    #         shape = 21, colour = "black", size = 2, stroke = 0.3, alpha=0.3) +
    geom_ribbon(aes(ymin = CI.lo_1, ymax = CI.hi_1),
                size=2, alpha = 0.1, linetype = 0) +
    #scale_color_manual(values=my_colours) +
    #scale_fill_manual(values=my_colours) +
    #scale_y_continuous(name=expression(~italic("Fv / Fm")),
    #                   limits = c(0.5, 0.61), 
    #                   breaks = seq(0.4, 0.6, by=0.1), expand = c(0,0))+
    #scale_x_continuous("Days in the experiment", limits = c(0, 78),
    #                 breaks = seq(0, 76, by=7), expand = c(0,0))+
    
     stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1,
                   position = position_dodge(1) )+
    stat_summary(fun.y=mean, geom="line", position = position_dodge(1), 
                  linetype=1, alpha=1) + 
    # stat_summary(fun.y=mean, geom="point", size =1,
    #             position=position_dodge(width=1), alpha=0.8)  +
    ggthe_bw
  
  Model_Ss_1b_plot + facet_wrap(Genotype~Community)

All phases, except baseline

  • Dasys as a factor
qPCR.variables_3<-qPCR.variables_2[(qPCR.variables_2$DaysF!="0"), ]  

Sy.Summary <- plyr::ddply (qPCR.variables_3, . (DaysF, Community, Treatment),
              summarise,
              Sy_mean = mean (TotalSH, na.rm = T),
              Sy_sd = sd (TotalSH, na.rm = T))
Sy.Summary
write.csv(Sy.Summary, "Outputs/meanSH.csv")
Treatment<- ggplot(
    qPCR.variables_3, aes(x = Days, y = logSH, 
                          colour = Treatment, shape = Treatment)) +
  Fill.colour+
    stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1,
                   position = position_dodge(3) )+
    stat_summary(fun.y=mean, geom="line", position = position_dodge(3), 
                  linetype=1, alpha=1)+
     stat_summary(fun.y=mean, geom="point", size =1,
                 position=position_dodge(3), alpha=0.8)  +
    ggthe_bw + facet_grid(~Community) +
   scale_y_continuous(breaks = seq(-5.5, 0.5, 0.5),
                    expand = c(0.1,0.1),
                  name=("log 10 (Total S/H cell ratio)")) +
  scale_x_continuous(name="Days in the experiment",
                         limits = c(75,113),
                         breaks = seq(75, 110, 15),  
                         expand = c(0, 0))  +
  theme(legend.position = c(.80, 0.24))+
    annotate("text", x = 78, y = -5, label = "(C)", size=3, colour="gray")+
    annotate("text", x = 100, y = -5, label = "(H)", size=3, colour="gray")
Treatment

#ggsave(file="Outputs/TotSH_Community.svg", plot=Treatment, width=4.5, height=3.5)
LM_Ssid.qpCR<-lmer(logSH ~ Treatment * DaysF * Community + (1|Replicate) + (1|Genotype/Fragment), 
                   data= qPCR.variables_3)
    step(LM_Ssid.qpCR)
## Backward reduced random-effect table:
## 
##                         Eliminated npar  logLik    AIC    LRT Df
## <none>                               31 -220.31 502.61          
## (1 | Fragment:Genotype)          1   30 -220.31 500.61  0.000  1
## (1 | Replicate)                  0   29 -237.04 532.09 33.477  1
## (1 | Genotype)                   0   29 -249.84 557.68 59.068  1
##                         Pr(>Chisq)    
## <none>                                
## (1 | Fragment:Genotype)      0.997    
## (1 | Replicate)          7.210e-09 ***
## (1 | Genotype)           1.523e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                           Eliminated  Sum Sq Mean Sq NumDF  DenDF F value
## Treatment:DaysF:Community          1  2.3236  0.2905     8 314.52  1.6881
## Treatment:Community                2  0.8550  0.2138     4 322.58  1.2212
## Treatment:DaysF                    0  3.2886  0.8221     4 326.73  4.6834
## DaysF:Community                    0 23.9600  5.9900     4 326.81 34.1226
##                              Pr(>F)    
## Treatment:DaysF:Community  0.100414    
## Treatment:Community        0.301634    
## Treatment:DaysF            0.001092 ** 
## DaysF:Community           < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## logSH ~ Treatment + DaysF + Community + (1 | Replicate) + (1 | 
##     Genotype) + Treatment:DaysF + DaysF:Community
    anova(LM_Ssid.qpCR)
    ranova(LM_Ssid.qpCR)
LM_Ssid.qpCR<-lmer(logSH ~ Treatment + DaysF + Community + 
                     (1 | Replicate) + (1 | Genotype) + 
                     Treatment:DaysF + DaysF:Community, data= qPCR.variables_3)
      step(LM_Ssid.qpCR)
## Backward reduced random-effect table:
## 
##                 Eliminated npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>                       18 -222.50 480.99                         
## (1 | Replicate)          0   17 -237.29 508.58 29.595  1  5.325e-08 ***
## (1 | Genotype)           0   17 -251.70 537.40 58.408  1  2.130e-14 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                 Eliminated  Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)
## Treatment:DaysF          0  3.2886  0.8221     4 326.73  4.6834  0.001092
## DaysF:Community          0 23.9600  5.9900     4 326.81 34.1226 < 2.2e-16
##                    
## Treatment:DaysF ** 
## DaysF:Community ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## logSH ~ Treatment + DaysF + Community + (1 | Replicate) + (1 | 
##     Genotype) + Treatment:DaysF + DaysF:Community
      anova(LM_Ssid.qpCR)
      ranova(LM_Ssid.qpCR)
      summary(LM_Ssid.qpCR)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: logSH ~ Treatment + DaysF + Community + (1 | Replicate) + (1 |  
##     Genotype) + Treatment:DaysF + DaysF:Community
##    Data: qPCR.variables_3
## 
## REML criterion at convergence: 445
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4207 -0.6488  0.1561  0.6873  2.2624 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  Genotype  (Intercept) 0.3836   0.6193  
##  Replicate (Intercept) 0.0364   0.1908  
##  Residual              0.1755   0.4190  
## Number of obs: 349, groups:  Genotype, 7; Replicate, 2
## 
## Fixed effects:
##                         Estimate Std. Error         df t value Pr(>|t|)
## (Intercept)            -1.646837   0.300750   8.678241  -5.476 0.000445
## TreatmentN              0.026826   0.083953 326.392615   0.320 0.749526
## TreatmentN+P           -0.155045   0.084638 326.495151  -1.832 0.067885
## DaysF100               -1.047853   0.137310 326.459574  -7.631 2.58e-13
## DaysF111               -1.308796   0.141498 326.622442  -9.250  < 2e-16
## CommunityC3             1.660629   0.206771 291.644381   8.031 2.40e-14
## CommunityD              0.199342   0.138081 332.129562   1.444 0.149778
## TreatmentN:DaysF100     0.394931   0.129108 326.496928   3.059 0.002405
## TreatmentN+P:DaysF100   0.313005   0.134932 326.499993   2.320 0.020972
## TreatmentN:DaysF111     0.003397   0.133201 327.432319   0.026 0.979667
## TreatmentN+P:DaysF111   0.332339   0.137490 326.610740   2.417 0.016188
## DaysF100:CommunityC3   -0.492060   0.152909 326.388468  -3.218 0.001420
## DaysF111:CommunityC3   -1.576303   0.176948 326.557888  -8.908  < 2e-16
## DaysF100:CommunityD    -0.126905   0.138844 326.531911  -0.914 0.361387
## DaysF111:CommunityD     0.140488   0.141114 327.593452   0.996 0.320197
##                          
## (Intercept)           ***
## TreatmentN               
## TreatmentN+P          .  
## DaysF100              ***
## DaysF111              ***
## CommunityC3           ***
## CommunityD               
## TreatmentN:DaysF100   ** 
## TreatmentN+P:DaysF100 *  
## TreatmentN:DaysF111      
## TreatmentN+P:DaysF111 *  
## DaysF100:CommunityC3  ** 
## DaysF111:CommunityC3  ***
## DaysF100:CommunityD      
## DaysF111:CommunityD      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Emmeans
Ssid.YII.emm<-emmeans(LM_Ssid.qpCR, ~Treatment*DaysF*Community)
     #contrast(Ssid.YII.emm, "tukey")

# Pairwise comparisons
    #pairs(Ssid.YII.emm) # same than contrast(Ssid.YII.emm, "tukey")
  plot(emmeans(LM_Ssid.qpCR, ~Treatment|DaysF), comparisons = TRUE) +
        theme_bw() # Tukey comparission (do not trust CI to compare EMMs)

    plot(emmeans(LM_Ssid.qpCR, ~Treatment|DaysF|Community), comparisons = TRUE) +
      coord_flip(xlim = NULL, ylim = NULL, expand = TRUE) + facet_grid(Community~DaysF) +
      theme_bw()

  #CLD  
  SH_Multicomp<-cld(Ssid.YII.emm, by=NULL) # compact-letter display
      SH_Multicomp<-SH_Multicomp[order(SH_Multicomp$DaysF, 
                                       SH_Multicomp$Community,
                                       SH_Multicomp$Treatment), ]
    #write.csv(SH_Multicomp, "Outputs/Ssid_SH_Multicomp.csv")  
 # 2. Predict values:
    pred_Ssid1 <- predict(LM_Ssid.qpCR, re.form = NA)

  #3. Bootstrap CI:
    Ss.boot1 <- bootMer(LM_Ssid.qpCR, predict, nsim = 1000, re.form = NULL) # include random effects, reduce CI lot!
    std.err <- apply(Ss.boot1$t, 2, sd)
    CI.lo_1 <- pred_Ssid1 - std.err*1.96
    CI.hi_1 <- pred_Ssid1 + std.err*1.96

  #Plot
  Model_Ss_1b_plot<- ggplot(
    qPCR.variables_3, aes(x = Days, y = logSH, colour = Treatment)) +
    geom_line(aes(y = pred_Ssid1),size=2) +
    #geom_point(aes(fill=factor(Treatment)),
    #         shape = 21, colour = "black", size = 2, stroke = 0.3, alpha=0.3) +
    geom_ribbon(aes(ymin = CI.lo_1, ymax = CI.hi_1),
                size=2, alpha = 0.1, linetype = 0) +
    #scale_color_manual(values=my_colours) +
    #scale_fill_manual(values=my_colours) +
    #scale_y_continuous(name=expression(~italic("Fv / Fm")),
    #                   limits = c(0.5, 0.61), 
    #                   breaks = seq(0.4, 0.6, by=0.1), expand = c(0,0))+
    #scale_x_continuous("Days in the experiment", limits = c(0, 78),
    #                 breaks = seq(0, 76, by=7), expand = c(0,0))+
    
     stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1,
                   position = position_dodge(1) )+
    stat_summary(fun.y=mean, geom="line", position = position_dodge(1), 
                  linetype=1, alpha=1)+
    # stat_summary(fun.y=mean, geom="point", size =1,
    #             position=position_dodge(width=1), alpha=0.8)  +
    ggthe_bw
  
  Model_Ss_1b_plot + facet_wrap(~Community)

  • Days as continous
LM_Ssid.qpCR2<-lmer(logSH ~ Treatment * Days * Community + (1|Replicate) + (1|Genotype/Fragment), 
                   data= qPCR.variables_2)
    step(LM_Ssid.qpCR2)
## Backward reduced random-effect table:
## 
##                         Eliminated npar  logLik    AIC   LRT Df Pr(>Chisq)
## <none>                               22 -568.24 1180.5                    
## (1 | Replicate)                  1   21 -568.24 1178.5  0.00  1          1
## (1 | Fragment:Genotype)          2   20 -568.24 1176.5  0.00  1          1
## (1 | Genotype)                   0   19 -589.20 1216.4 41.93  1  9.461e-11
##                            
## <none>                     
## (1 | Replicate)            
## (1 | Fragment:Genotype)    
## (1 | Genotype)          ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                          Eliminated  Sum Sq Mean Sq NumDF  DenDF F value
## Treatment:Days:Community          1  1.7016  0.4254     4 461.15  0.9054
## Days:Community                    2  1.7500  0.8750     2 465.80  1.8642
## Treatment:Community               3  3.2309  0.8077     4 467.27  1.7133
## Community                         0 22.9355 11.4678     2 364.82 24.1860
## Treatment:Days                    0  3.0544  1.5272     2 471.54  3.2209
##                             Pr(>F)    
## Treatment:Days:Community    0.4606    
## Days:Community              0.1562    
## Treatment:Community         0.1458    
## Community                1.369e-10 ***
## Treatment:Days              0.0408 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## logSH ~ Treatment + Days + Community + (1 | Genotype) + Treatment:Days
    anova(LM_Ssid.qpCR2)
    ranova(LM_Ssid.qpCR2)
LM_Ssid.qpCR2<-lmer(logSH ~Treatment + Days + Community + (1|Genotype) + Treatment:Days,
                   data= qPCR.variables_2)
      step(LM_Ssid.qpCR2)
## Backward reduced random-effect table:
## 
##                Eliminated npar  logLik    AIC    LRT Df Pr(>Chisq)    
## <none>                      10 -542.81 1105.6                         
## (1 | Genotype)          0    9 -564.04 1146.1 42.471  1  7.173e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Backward reduced fixed-effect table:
## Degrees of freedom method: Satterthwaite 
## 
##                Eliminated  Sum Sq Mean Sq NumDF  DenDF F value    Pr(>F)
## Community               0 22.9355 11.4678     2 364.82 24.1860 1.369e-10
## Treatment:Days          0  3.0544  1.5272     2 471.54  3.2209    0.0408
##                   
## Community      ***
## Treatment:Days *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Model found:
## logSH ~ Treatment + Days + Community + (1 | Genotype) + Treatment:Days
      anova(LM_Ssid.qpCR2)
      ranova(LM_Ssid.qpCR2)
      summary(LM_Ssid.qpCR2)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## logSH ~ Treatment + Days + Community + (1 | Genotype) + Treatment:Days
##    Data: qPCR.variables_2
## 
## REML criterion at convergence: 1085.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.5286 -0.6117 -0.1087  0.6507  2.3774 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  Genotype (Intercept) 0.6803   0.8248  
##  Residual             0.4741   0.6886  
## Number of obs: 486, groups:  Genotype, 7
## 
## Fixed effects:
##                     Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)        -1.853179   0.360321   8.344461  -5.143 0.000773 ***
## TreatmentN         -0.096926   0.136250 471.094970  -0.711 0.477201    
## TreatmentN+P       -0.266358   0.143608 471.050405  -1.855 0.064257 .  
## Days               -0.008168   0.001209 471.227691  -6.757 4.18e-11 ***
## CommunityC3         1.818818   0.283439 290.421477   6.417 5.64e-10 ***
## CommunityD          0.262047   0.169203 453.518949   1.549 0.122148    
## TreatmentN:Days     0.003229   0.001701 471.924639   1.898 0.058301 .  
## TreatmentN+P:Days   0.004293   0.001799 471.108284   2.386 0.017420 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmnN TrtN+P Days   CmmnC3 CmmntD TrtN:D
## TreatmentN  -0.183                                          
## TreatmntN+P -0.185  0.480                                   
## Days        -0.248  0.591  0.564                            
## CommunityC3 -0.371 -0.026  0.005  0.053                     
## CommunityD  -0.382 -0.009  0.011  0.047  0.581              
## TrtmntN:Dys  0.149 -0.834 -0.400 -0.706  0.063 -0.018       
## TrtmntN+P:D  0.141 -0.399 -0.839 -0.669  0.021  0.020  0.477
 # 2. Predict values:
    pred_Ssid1 <- predict(LM_Ssid.qpCR2,re.form = NA)

  #3. Bootstrap CI:
    Ss.boot1 <- bootMer(LM_Ssid.qpCR2, predict, nsim = 1000, re.form = NULL) # include random effects, reduce CI lot!
    std.err <- apply(Ss.boot1$t, 2, sd)
    CI.lo_1 <- pred_Ssid1 - std.err*1.96
    CI.hi_1 <- pred_Ssid1 + std.err*1.96

  #Plot
  Model_Ss_1b_plot<- ggplot(
    qPCR.variables_2, aes(x = Days, y = logSH, colour = Treatment)) +
    geom_line(aes(y = pred_Ssid1),size=2) +
    #geom_point(aes(fill=factor(Treatment)),
    #         shape = 21, colour = "black", size = 2, stroke = 0.3, alpha=0.3) +
    geom_ribbon(aes(ymin = CI.lo_1, ymax = CI.hi_1),
                size=2, alpha = 0.1, linetype = 0) +
    #scale_color_manual(values=my_colours) +
    #scale_fill_manual(values=my_colours) +
    #scale_y_continuous(name=expression(~italic("Fv / Fm")),
    #                   limits = c(0.5, 0.61), 
    #                   breaks = seq(0.4, 0.6, by=0.1), expand = c(0,0))+
    #scale_x_continuous("Days in the experiment", limits = c(0, 78),
    #                 breaks = seq(0, 76, by=7), expand = c(0,0))+
    
     stat_summary(fun.data = "mean_cl_boot",geom = "errorbar", width = 1,
                   position = position_dodge(1) )+
    stat_summary(fun.y=mean, geom="line", position = position_dodge(1), 
                  linetype=1, alpha=1) + 
    # stat_summary(fun.y=mean, geom="point", size =1,
    #             position=position_dodge(width=1), alpha=0.8)  +
    ggthe_bw
  
  Model_Ss_1b_plot + facet_wrap(Genotype~Community)